rm(list=ls())
if(!require('pacman'))install.packages('pacman')
pacman::p_load(foreign,readstat13,dplyr,foreign,nnet, ggplot2, tidyr, stargazer, car, gmodels, Hmisc, pscl,
               scales, data.table, zoo, mice, MALDIquant, arsenal, countrycode, readxl, clusterSEs, rms, lfe, 
               ggeffects, sandwich, miceadds, here)

library(foreign)
library(readstata13)
library(dplyr)
library(foreign)
library(nnet)
library(ggplot2)
library(tidyr)
library(stargazer)
library(car)
library(gmodels)
library(Hmisc)
library(pscl)
library("scales")
library(data.table)
library(zoo)
library(mice)
library("MALDIquant")
library(arsenal)
library(countrycode)
library(readxl)
library(clusterSEs)
library(rms)
library(lfe)
library(ggeffects)
library("sandwich")
library("miceadds")
library(here)

setwd(this.path::here())

#####Load Data#####

Data <-read.csv("Replication_II2025/TIES_IdealPointMatched_TradeSanctions.csv") #For main analysis
Data_realizedsanction <- read.csv("Replication_II2025/RealizedSanctions_Formatted.csv") #Data for analysis of realized sanctions
load("Replication_II2025/HAGE.RData")  #Data for robustness test
Hagechance <- x #Rename data
TIES <- read.csv("Replication_II2025/TIESPrelim.csv")




####Descriptive Figures####

####Figure 1####
####Count Multi Sanctions####
TIES$multi <- ifelse(is.na(TIES$sender2), 0, 1)
count_multi <- TIES %>% group_by(startyear) %>% count(multi)
total_multi <- count_multi %>% group_by(startyear) %>% summarise(total = sum(n))
count_multi <- left_join(count_multi, total_multi, by = ("startyear"))
count_multi$percent <- count_multi$n/count_multi$total
count_multi_sub <- subset(count_multi, multi==1)
plot(count_multi_sub$startyear, count_multi_sub$percent, type = "h", xlab = "Year", ylab = "Percent of Multilateran Sanctions")


###Figure 3 -- Plotting Matched Countries###
##Ideal Point Distance Plot (Figure 3a)

plotdata <- cbind(Data$primarysender, Data$matchedcountryp,
                  Data$diffabspnorm) #Extract primary sender-matched country pairing
plotdata <- as.data.frame(plotdata) #Convert to dataframe 
colnames(plotdata) <- c("primarysender", "matchedcountry", "idealpointdiff") #Rename columns
plotdata <- plotdata %>% filter(!is.na(idealpointdiff)) #Drop NA's (these don't appear in the empirical analysis, either)


Fig3a <- plotdata %>%
  ggplot( aes(x=idealpointdiff)) + ggtitle("Distribution of Ideal Point Difference") + 
  xlab("Ideal point difference") + ylab("Count") +
  geom_density(fill="blue", color="#e9ecef", alpha=0.8) ###This creates Figure 3a

##RMP Added Mathced States Distribution (Figure 3b)##

Fig3b <- Data %>%
  ggplot( aes(x=relative_market_power_added)) + ggtitle("Distribution of Added Relative Market Power (RMP)") + 
  xlab("RMP added") + ylab("Density") +
  geom_density(fill="blue", color="#e9ecef", alpha=0.8) #This creates Figure 3b




####Analysis####
####Table 1####
fit1 <- glm(multi.x ~ diffabspnorm + relative_wealth + cinc+ target_export_portfolio , 
            data = Data, family = "binomial") 
summary(fit1)
a <- coeftest(fit1, vcov. = vcovCL(fit1, cluster = Data$primarysender, type = "HC0"))
summary(a)
stargazer(a) #Create Table 1

####Figure 4####
mydf <- ggpredict(fit1, terms = "diffabspnorm[all]")
Fig4 <- ggplot(mydf, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
  xlab("Policy Preference Distance") + ylab("Predicted Probability of Multilateral Sanction") + theme_bw() #Create Figure 4




####Table 3####
fit2 <- glm(multi.x~relative_market_power_added+ relative_wealth + target_export_portfolio  + cinc + factor(democ_bi), 
            data = Data, 
            family = "binomial")
summary(fit2)

fit2d <- glm(multi.x~relative_market_power_added+ relative_wealth + target_export_portfolio  + cinc , 
             data = Data, 
             family = "binomial")
summary(fit2d)

b <- coeftest(fit2, vcov. = vcovCL(fit2, cluster = Data$primarysender, type = "HC0"))
bd <- coeftest(fit2d, vcov. = vcovCL(fit2d, cluster = Data$primarysender, type = "HC0"))
b
bd
stargazer(b, bd) #This creates Table 3



####Figure 5####

mydf2 <- ggpredict(fit2, terms = "relative_market_power_added[all]")
Fig5 <- ggplot(mydf2, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
  xlab("Added Coercive Power") + ylab("Predicted Probability of Multilateral Sanction") + theme_bw()




####Table 4 -- Democracies as Primary Sender####
Data_democ <- subset(Data, democ_bi == 1) #Subset Data to only democratic primary senders 

fitl3c3 <- glm(multi.x~diffabspnorm + relative_wealth + cinc + target_export_portfolio , data = Data_democ, family = "binomial") #regression
summary(fitl3c3)
d <- coeftest(fitl3c3, vcov. = vcovCL(fitl3c3, cluster = Data_democ$primarysender, type = "HC0")) #cluster resudiual
d
stargazer(fitl3c3) #This creates Table 4



####Figure 6 -- Trade Objective vs. Non-Trade Objective Divergent Prferences####

Data$tradeissue <- as.factor(Data$tradeissue) #Code binary indicator for whether sanction objective is trade-related 
fit15 <- glm(multi.x ~ diffabspnorm + factor(tradeissue) + relative_wealth + cinc+ target_export_portfolio, data = Data, family = "binomial") #regression analysis
summary(fit15)


mydf1 <- ggeffects::ggpredict(fit15, terms = c("diffabspnorm[all]", "tradeissue")) #setting up data for graph creation
mydf1$group <- ifelse(mydf1$group == "0", "Non-Trade", "Trade")
names(mydf1)[names(mydf1) == "group"] <- "Objective"
aB1<- ggplot(mydf1, aes(x, predicted, color = Objective)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Objective), alpha = .1) +  geom_line(aes(colour = Objective),
                                                                                                 size = 1) +
  xlab("Policy Preference Distance") + ylab("Predicted Probability of Multilateral Sanction") + theme_bw() #Create Figure 6



####Figure 7 -- Trade vs Non-Trade Objective RMP####<- NEED TO RECEHCK THIS

fit55 <- glm(multi.x ~ relative_market_power_added + factor(tradeissue) + relative_wealth + cinc+ target_export_portfolio + factor(democ_bi), data = Data, family = "binomial") #regression
summary(fit55)

mydf3 <- ggpredict(fit55, terms = c("relative_market_power_added[all]", "tradeissue"))
mydf3$group <- ifelse(mydf3$group == "0", "Non-Trade", "Trade")
names(mydf3)[names(mydf3) == "group"] <- "Objective"
a2<- ggplot(mydf3, aes(x, predicted, color = Objective)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Objective), alpha = .1) +  geom_line(aes(colour = Objective),
                                                                                                 size = 1) +
  xlab("Added Coercive Leverage") + ylab("Predicted Probability of Multilateral Sanction") + theme_bw() 

a2


####Table 5 -- Realized Sanctions####

fit3 <- glm(chosen ~ relative_market_power+diffidnorm + cinc + cincr + shared_democ  , data = Data_realizedsanction, family = "binomial") 
summary(fit3)
yy <- coeftest(fit3, vcov. = vcovCL(fit3, cluster = Data_realizedsanction$caseid, type = "HC0"))
yy


fit31 <- glm(chosen ~ relative_market_power+diffidnorm  , data = Data_realizedsanction, family = "binomial") 
summary(fit31)
yy1 <- coeftest(fit31, vcov. = vcovCL(fit31, cluster = Data_realizedsanction$caseid, type = "HC0"))
yy1

fit32 <- glm(chosen ~ relative_market_power+diffidnorm +cinc , data = Data_realizedsanction, family = "binomial") 
summary(fit32)
yy2 <- coeftest(fit32, vcov. = vcovCL(fit32, cluster = Data_realizedsanction$caseid, type = "HC0"))
yy2

fit33 <- glm(chosen ~ relative_market_power+diffidnorm +shared_democ , data = Data_realizedsanction, family = "binomial") 
summary(fit33)
yy3 <- coeftest(fit33, vcov. = vcovCL(fit33, cluster = Data_realizedsanction$caseid, type = "HC0"))
yy3

fit34 <- glm(chosen ~ relative_market_power+diffidnorm +cincr , data = Data_realizedsanction, family = "binomial") 
summary(fit34)
yy4 <- coeftest(fit34, vcov. = vcovCL(fit34, cluster = Data_realizedsanction$caseid, type = "HC0"))
yy4



stargazer(yy1, yy2, yy3, yy4, yy) #This is used to create Table 5



####Robustness Check -- Table 2####
###Robustness using Chance-Corrected Cohn's Kappa###
#Get the needed variables 
Kappa <- cbind(Hagechance$year, Hagechance$ccode1, Hagechance$ccode2, Hagechance$kappavv)
Kappa <- as.data.frame(Kappa)
colnames(Kappa) <- c("year", "primarysender", "ccode2", "kappavv")
matchingprim <- cbind(Data$caseid, Data$year, Data$primarysender)
matchingprim <- as.data.frame(matchingprim)
colnames(matchingprim) <- c("caseid", "year", "primarysender")

#remove rows with duplicate observations
Kappa <- Kappa %>% 
  as_tibble() %>% 
  mutate(duplicates = if_else(primarysender == ccode2,
                              TRUE,
                              FALSE)) %>% 
  filter(duplicates == FALSE)

Kappa$duplicates <- NULL
Kappa <- filter(Kappa, year > 1960)


matchingprim <- left_join(matchingprim, Kappa, by = c("year", "primarysender"))
matchingprim <- filter(matchingprim, year > 1960)

testing <- matchingprim %>%
  group_by(caseid, year, primarysender) %>%
  slice(which.min(abs(kappavv - 1)))

colnames(testing) <- c("caseid", "year", "primarysender", "matched_hage", "kappavv")


Data <- left_join(Data, testing, by = c("year", "primarysender", "caseid"))



#Testing

fith <- glm(multi.x ~ kappavv + relative_wealth + cinc+ target_export_portfolio , 
            data = Data, family = "binomial")
summary(fith)
stargazer(fith)


fithb <- coeftest(fith, vcov. = vcovCL(fith, cluster = Data$primarysender, type = "HC0"))
stargazer(fithb)
fithb


####Robustness using Chance-Corrected Scott's Pi####

Pi <- cbind(Hagechance$year, Hagechance$ccode1, Hagechance$ccode2, Hagechance$pivv)
Pi <- as.data.frame(Pi)
colnames(Pi) <- c("year", "primarysender", "ccode2", "pivv")
matchingprim2 <- cbind(Data$caseid, Data$year, Data$primarysender)
matchingprim2 <- as.data.frame(matchingprim2)
colnames(matchingprim2) <- c("caseid", "year", "primarysender")

#remove rows with duplicate observations
Pi <- Pi %>% 
  as_tibble() %>% 
  mutate(duplicates = if_else(primarysender == ccode2,
                              TRUE,
                              FALSE)) %>% 
  filter(duplicates == FALSE)

Pi$duplicates <- NULL
Pi <- filter(Pi, year > 1960)


matchingprim2 <- left_join(matchingprim2, Pi, by = c("year", "primarysender"))
matchingprim2 <- filter(matchingprim2, year > 1960)

testingpi <- matchingprim2 %>%
  group_by(caseid, year, primarysender) %>%
  slice(which.min(abs(pivv - 1)))

colnames(testingpi) <- c("caseid", "year", "primarysender", "matched_hage_spi", "pivv")


Data <- left_join(Data, testingpi, by = c("year", "primarysender", "caseid"))

#Testing

fitpi <- glm(multi.x ~ pivv + relative_wealth + cinc+ target_export_portfolio , 
             data = Data, family = "binomial")
summary(fitpi)
stargazer(fitpi)


fithbpi <- coeftest(fitpi, vcov. = vcovCL(fitpi, cluster = Data$primarysender, type = "HC0"))
stargazer(fithb, fithbpi) ###This creates Table 2








